install.packages("pacman")
pacman::p_load(tidyverse, survival, readxl, ggplot2, survminer)
data <- read_excel("data.xlsx")
# Remove NA values
data$export_city <- as.numeric(data$export_city)
data$disaster <- as.numeric(data$disaster)
data <- na.omit(data[, c("city","population", "democratic", "density", "export_city", "disaster",
"risk", "readiness", "coastal", "t2", "gcom", "airport")])
# Suppressing unnecessary messages
options(warn=-1)
# Declare survival dataset
data$surv <- ifelse(data$gcom == 1 & data$t2 < 11, 1, 0)
# descriptive statistics
summary(data)
summary(data$export_city) / 1000000
data$population <- log(data$population + 1)
data$density <- log(data$density + 1)
data$export_city <- log(data$export_city + 1)
##Analysis
# Fit the Cox model
cox_model1 <- coxph(Surv(t2, surv) ~ readiness, data = data)
cox_model2 <- coxph(Surv(t2, surv) ~ readiness + population  + density + democratic +
coastal + risk + export_city, data = data)
cox_model3 <- coxph(Surv(t2, surv) ~ readiness + population  + density + democratic +
coastal + risk + export_city + frailty(city), data = data)
# PH assumption
PH <- cox.zph(cox_model2, transform="km")
PH
# Extract the summary of the Cox model
cox_summary1 <- summary(cox_model1)
cox_summary2 <- summary(cox_model2)
cox_summary3 <- summary(cox_model3)
cox_summary1
cox_summary2
cox_summary3
extractAIC(cox_model1)
extractAIC(cox_model2)
extractAIC(cox_model3)
## counterfactual scenarios
library(MASS); library(simcf); library(tile)
# for simcf and tile, please refer to: https://faculty.washington.edu/cadolph/index.php?page=60
simCoxPH <- function(sims, before, after, baseline=0) {
crossRR <- sims^(after - before)
res <- c(mean(crossRR), quantile(crossRR, probs=c(0.025, 0.975)))
res
}
sims <- 10000
simbetas <- exp(mvrnorm(sims, coef(cox_model2), vcov(cox_model2)))
CFnames <- CFs <- NULL
CFnames <- c(CFnames, "Scenario 1")
CFs <- rbind(CFs,
simCoxPH(simbetas[,1],
mean(data$readiness),
mean(data$readiness) + sd(data$readiness)
)
)
CFnames <- c(CFnames, "Scenario 2")
CFs <- rbind(CFs,
simCoxPH(simbetas[,2],
mean(data$population),
mean(data$population) + 1
)
)
CFnames <- c(CFnames, "Scenario 3")
CFs <- rbind(CFs,
simCoxPH(simbetas[,4],
0,
1
)
)
CFnames <- c(CFnames, "Scenario 4")
CFs <- rbind(CFs,
simCoxPH(simbetas[,7],
mean(data$export_city),
mean(data$export_city) + 1
)
)
trace <- ropeladder(x=CFs[,1],
lower=CFs[,2],
upper=CFs[,3],
labels=CFnames,
entryheight=0.2,
col=c("#e6ab2b", "#8ee62b", "#2b83e6", "#aa2be6"),
size=0.5, lex=1.5,
lineend="square",
plot=1
)
vertmark <- linesTile(x=c(1,1), y=c(0,1), layer=12, plot=1)
file <- "Figure1"
tile(trace, vertmark,
output=list(outfile=file, width=7.0),
limits=c(0, 20),
gridlines=list(type="xt"),
topaxis=list(add=T, log=F, at=c(1, 5, 10, 15, 20),
labels=c("1x", "5x", "10x", "15x", "20x")),
xaxis=list(log=F, at=c(1, 5, 10, 15, 20),
labels=c("1x", "5x", "10x", "15x", "20x")),
width=list(plot=2),
height=list(topaxistitle=2)
)
tile(trace, vertmark,
limits=c(0, 20),
gridlines=list(type="xt"),
topaxis=list(add=TRUE, log=TRUE, at=c(0, 1, 2, 5, 10, 20),
labels=c("0x", "1x", "10x", "15x", "20x")),
xaxis=list(log=TRUE, at=c(0, 5, 10, 15, 20),
labels=c("0x", "5x", "10x", "15x", "20x")),
topaxistitle=list(labels="Relative risk of joining GCoM", y=0.8),
width=list(plot=2),
height=list(topaxistitle=2),
output=list(outfile=file,width=7.0)
) # change the font later
